class BlockEP_t
{
	public:
		BlockEP_t(mat L, mat b)
		{
			_L=L;
			_p=L.n_cols;
			d=new boost::math::normal_distribution<>(0,1); 
			boost::random::uniform_01<> U;
			_runif=new RandomG::Random<boost::random::uniform_01<> >(U);
			_b=b;
			mat t(2,1);
			G=t;
			mat I(2,2);
			I.fill(1);
			_I=diagmat(I);
		}
		void operator()(int j, mat x,int k, double *z, double *w)
		{
			if(j<_p-1)
			{
				double tr1=_b(0,j);
				double tr2=_b(0,j+1);
				for(int i=0;i<j;i++) 
				{
					tr1-=_L(j,i)*x(0,i);
					tr2-=_L(j+1,i)*x(0,i);
				}
				tr1=tr1/_L(j,j);
				if(k==0)
				{
					G(0,0)=_L(j+1,j);
					G(1,0)=_L(j+1,j+1);
					norm=sqrt(dot(G,G));
					norm3=norm*norm*norm;
				}
				double cd= cdf(*d,tr2/norm);
				if(cd > .999999999) cd=.999999999;
				if(cd < .0000000001) cd=.0000000001;
				mat dlZm=((double)1/(1-cd)*pdf(*d,tr2/norm))*G/norm;
				mat dlZS=((double)1/(1-cd)*pdf(*d,tr2/norm)*tr2/(norm3))*G*G.t();


				mat Sig=_I-dlZm*dlZm.t()+2*dlZS;
				double sig=sqrt(Sig(0,0));
				double m=dlZm(0,0);	
				
				
				double cd2=cdf(*d,tr1);	
				double cd4=cdf(*d,(tr1-m)/sig);	
				double u=(*_runif)(); 
				double arg=u+(1.-u)*cd4;
				if(arg > .999999999) arg=.999999999;
				if(arg < .0000000001) arg=.0000000001;
				*z = quantile(*d,arg)*sig+m;
				double tpz=(tr2-_L(j+1,j)*(*z))/_L(j+1,j+1);
				
				double cd3=cdf(*d,tpz);	
				*w=(1-cd3)*pdf(*d,*z)*sig*(1-cd4)/(pdf(*d,(*z-m)/sig)*(1-cd2));
				//cout << *w << " ";
			}else{
				double mu=0.0;for(int kk=0;kk<j;kk++) mu += _L(j,kk)*x(0,kk);
				double tpz = (_b(0,j)-mu)/_L(j,j);
				
				//cout << "tpz " << tpz;
				double pa = cdf(*d,tpz);
				//cout << prod[i] << " ";
				double u =(*_runif)(); 
				//u=Nt[kk];
				double arg=u+(1.-u)*pa;
				if(arg > .999999999) arg=.999999999;
				if(arg < .0000000001) arg=.0000000001;
				*w=1;
				*z = quantile(*d,arg);
			}

		}
	private:
		mat _L;
		int _p;
		mat _b;
		mat _I;
		boost::math::normal_distribution<> *d;
		RandomG::Random<boost::random::uniform_01<> > *_runif;
		mat G;
		double norm3;
		double norm;
};
